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Abstract. 

In this paper, we address the accuracy of the results for the overdetermined full rank linear least squares 
problem. We recall theoretical results obtained in [2] on conditioning of the least squares solution and the 
components of the solution when the matrix perturbations are measured in Frobenius or spectral norms. Then we 
define computable estimates for these condition numbers and we interpret them in terms of statistical quantities. 
In particular, we show that, in the classical linear statistical model, the ratio of the variance of one component of 
the solution by the variance of the right-hand side is exactly the condition number of this solution component when 
perturbations on the right-hand side are considered. We also provide fragment codes using LAPACK ^ routines 
to compute the variance-covariance matrix and the least squares conditioning and we give the corresponding 
computational cost. Finally we present a small historical numerical example that was used by Laplace |19j for 
computing the mass of Jupiter and experiments from the space industry with real physical data. 
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1. Introduction. We consider the linear least squares problem (LLSP) mina-gRn \\Ax — &II2, 
where b £ R™ and A G ]U™x" is a matrix of full column rank n. 

Our concern comes from the following observation: in many parameter estimation problems, there 
may be random errors in the observation vector b due to instrumental measurements as well as 
roundoff errors in the algorithms. The matrix A may be subject to errors in its computation 
(approximation and/or roimdoff errors). In such cases, while the condition number of the matrix 
A provides some information about the sensitivity of the LLSP to perturbations, a single global 
conditioning quantity is often not relevant enough since we may have significant disparity between 
the errors in the solution components. We refer to the last section of the manuscript for illustrative 
examples. 

There are several results for analyzing the accuracy of the LLSP by components. For linear systems 
Ax = b and for LLSP, [7] defines so called componentwise condition numbers that correspond to 
amplification factors of the relative errors in solution components due to perturbations in data 
A or b and explains how to estimate them. For LLSP, |17j proposes to estimate componentwise 
condition numbers by a statistical method. More recently, [2] developed theoretical results on 
conditioning of linear functionals of LLSP solutions. 

The main objective of our paper is to provide computable quantities of these theoretical values in 
order to assess the accuracy of an LLSP solution or some of its components. To achieve this goal, 
traditional tools for the numerical linear algebra practitioner are condition numbers or backward 
errors whereas the statistician usually refers to variance or covariance. Our purpose here is to show 
that these mathematical quantities coming either from numerical analysis or statistics are closely 
related. In particular, we will show in Equation (|3.3|1 that, in the classical linear statistical model, 
the ratio of the variance of one component of the solution by the variance of the right-hand side 
is exactly the condition number of this component when perturbations on the right-hand side are 
considered. In that sense, we attempt to clarify, similarly to |15| . the analogy between quantities 
handled by the linear algebra and the statistical approaches in linear least squares. Then we 
define computable estimates for these quantities and explain how they can be computed using the 
standard hbraries LAPACK or ScaLAPACK. 
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This paper is organized as follows. In Section [21 we recall and exploit some results of practical 
interest coming from [2]. We also define the condition numbers of an LLSP solution or one 
component of it. In Section [31 we recall some definitions and results related to the linear statistical 
model for LLSP, and we interpret the condition numbers in terms of statistical quantities. In 
Section [H we provide practical formulas and FORTRAN code fragments for computing the 
variance-covariance matrix and LLSP condition numbers using LAPACK. In Section[Sl we propose 
two numerical examples that show the relevance of the proposed quantities and their practical 
computation. The first test case is a historical example from Laplace and the second exam- 
ple is related to gravity field computations. Finally some concluding remarks are given in Section[Sl 

Throughout this paper we will use the following notations. We use the Frobenius norm ||.||^ 
and the spectral norm ||.||2 on matrices and the usual Euclidean norm jj.jjj on vectors. denotes 
the Moore-Penrose pseudo inverse of A, the matrix / is the identity matrix and is the i-th 
canonical vector. 

2. Theoretical background for linear least squares conditioning. Following the nota- 
tions in [2], we consider the function 

g : M™^" X M™ — > R'= , , 

A,b ^ giA,b) = L'^x{A,b)^L^{A'^A)-^A'^b, ^^■^> 

where L is an n x fc matrix, with k < n. Since A has full rank n, g is continuously F-differentiable 
in a neighbourhood of (A, b) and we denote by g' its F-derivative. 

Let a and /3 be two positive real numbers. In the present paper we consider the Euclidean norm 
for the solution space M*^. For the data space K™^" x R"*, we use the product norms defined by 



\\{A,b)hor2^^Ja^A\\l^^,+l3^\\b\\l a,/?>0. 

Following [lOj . the absolute condition number of g at the point {A,b) using the product norm 
defined above is given by: 

\W{A,b).{AA,Ab)\\, 

Kg,F or 2{A, b) = max 



(AA,Ab) |j(AA, A6)||f or 2 

The corresponding relative condition number of g at [A, b) is expressed by 



(rel) ( A u\ _ l^g.pjA, b) \\ (A, 6)||f pr 2 
or 2i^^f>) - ^{1^1 ■ 
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To address the special cases where only A (resp. b) is perturbed, we also define the quantities 

, ^, |||4.(A,b).AA|| , ,,, |j^(A,6).A&|| , 

Kg,F or 2iA) = maxAA ||AA|iF J ^"^^^"P' '^saW = maxAb ||Ab||, )■ 

Remark 1. The product norm for the data space is very flexible; the coefficients a and /3 
allow us to monitor the perturbations on A and b. For instance, large values of a (resp. /? ) enable 
us to obtain condition number problems where mainly b (resp. A) are perturbed. In particular, 
we will address the special cases where only b (resp. A) is perturbed by choosing the a and (3 
parameters as a = +00 and /? = 1 (resp. a = 1 and (3 = +00) since we have 

lim Kg.F or 2{A, b) = \ng,2{b) and lim k^.f or 2{A, b) = -%.f or 2{A). 

Q— »+oo ■ p ' p^+oo a 
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This can be justified as follows: 



Kg,F or 2{A b) = max 

{AA,Ab) 



max 

(AA.Ab) 



|f(A,6).A^ + |f(A,6).A6 



a2||AA|||„,2+/32||A6||^ 



1AA|1|„,, + |1A6| 



The above expression represents the operator norm of a linear functional depending continuously 
on a, and then we get 



lim Kg F or 2{-A,b) = max 

Q^+oo ' (AA.Ab) ^ /||^^||2 



F or 2 



|Afe| 



= = max 

2 Af) 



^i^gAb)- 



The proof is the same for the case where /? = +00. 



The condition numbers related to L'^xlA, b) are referred to as partial condition num- 
bers (PCN) of the LLSP with respect to the linear operator L in [2]. 

We are interested in computing the PCN for two special cases. The first case is when L is the 
identity matrix (conditioning of the solution) and the second case is when L is a canonical vector 
Ci (conditioning of a solution component). We can extract from [2] two theorems that can lead to 
computable quantities in these two special cases. 

Theorem 1. In the general case where {L € R"^*^), the absolute condition numbers of 
g{A, b) = L'^x{A, b) in the Frobenius and spectral norms can be respectively bounded as follows 

-^f{A,b)<KgAA,b)<f{A,b) 



where 



-^f{A,b)<Kg,2{A,b)<V2f{A,b) 



fiA,b)^ \\L^iA^A)-H'M 




(2.2) 



Theorem 2. In the two particular cases: 

1. L is a vector (L € R"j, or 

2. L is the n-by-n identity matrix (L — I) 

the absolute condition number of g{A, b) — L'^x{A, b) in the Frobenius norm is given by the formula: 

-.MA, b) = (\\L-iA-Ar\\l ^ + WL-At, (M + J_)y . 

Theorem [2] provides the exact value for the condition number in the Frobenius norm for our two 
cases of interest {L — Ci and L ~ I). From Theorem [1] we observe that 

4=Ks,f(A, b) < Kg,2(A, b) < V6Kg,F{A, b). (2.3) 
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which states that the partial condition number in spectral norm is of the same order of magnitude 
as the one in Frobenius norm. In the remainder of the paper, the focus is given to the partial 
condition number in Frobenius norm only. 

For the case L = I, the result of Theorem [2] is similar to [TT] and [TOl p. 92]. The upper bound 
for K2,F{A,b) that can be derived from Equation (|2.3p is also the one obtained by [TU] when we 
consider pertubations in A. 

Let us denote by b) the condition number related to the component x,; in Frobenius norm (i.e 
Ki{A, b) = Kg^p{A, b) where g{A^ b) — efx{A, b) = Xi{A, b)). Then replacing L by in Theorem[5] 
provides us with an exact expression for computing K.i{A, b), this gives 



.,:(A6)= lief + lief At||(^ + -^) . (2.4) 




Ki{A,b) will be referred to as the condition number of the solution component Xi. 

Let us denote by n^si^y b) the condition number related to the solution x in Frobenius norm 
(i.e KLs{A,b) = Kg^p{A,b) where g{A,b) — x{A,b)). Then Theorem [5] provides us with an exact 
expression for computing KLs{^,b), that is 




KLs{A,b)= UA'Ar^ r^ [11^ +-5^ . (2.5) 
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where we have used the fact that || (A"^^)"-^ II2 = ||^^||2- 
KLs(^:^) will be referred to as the condition number of the least squares solution. 

Note that [5] defines condition numbers for both x and r in order to derive error bounds for x and 
r but uses infinity-norm to measure perturbations. 

In this paper, we will also be interested in the special case where only b is perturbed (a = +00 
and (3 = 1). In this case, we will call Ki{b) the condition number of the solution component Xi, and 
KLs(^) the condition number of the least squares solution. When we restrict the perturbations to 
be on b, Equation (j2.4p simplifies to 



K,ib) = \\efA^\^, (2.6) 

and Equation ()2.5p simplifies to 

^Ls{b) = \\A^\.^. (2.7) 
This latter formula is standard and is in accordance with p. 29]. 
3. Condition numbers and statistical quantities. 

3.1. Background for the linear statistical model. We consider here the classical linear 
statistical model 

b = Ax + e, Ae R™^", 6 e M™, rank(A) = n, 

where e is a vector of random errors having expected value E(e) = and variance-covariance 
V{e) = afl. In statistical language, the matrix A is referred to as the regression matrix and the 
unknown vector x is called the vector of regression coefficients. 

Following the Gauss-Markov theorem [50], the least squares estimates x is the linear unbiased 
estimator of x satisfying 

\\Ax — 6II2 = min \\Ax — felU, 

with minimum variance-covariance equal to 

C = aUA^A)-\ (3.1) 



Moreover ||&— Ax\\2 is an unbiased estimate of cr^. This quantity is sometimes called the 
mean squared error (MSE). 

The diagonal elements cu of C give the variance of each component Xi of the solution. The off- 
diagonal elements Cy , i ^ j give the covariance between Xi and Xj . 
We define ax. as the standard deviation of the solution component Xi and we have 

CTxi = \fCii- (3.2) 

In the next section, we will prove that the condition numbers h) and ^^^(A, 6) can be related 
to the statistical quantities erf. and Uh. 

3.2. Perturbation on h only. Using Equation p.ip . the variance ca of the solution compo- 
nent Xi can be expressed as 

We note that {A^A)-^ = Atyl^^ so that 

Cu=alef{A^A^^)e, = al\\efA^\l. 

Using Equation p.2p . we get 
Finally from Equation (|2.6p . we get 

c^xi = crbHi(h). (3.3) 

Equation (|3.3p shows that the condition number (6) relates linearly the standard deviation of ct^ 
with the standard deviation of cr^ . . 

Now if we consider the constant vector £ of size n, we have (see [2U] ) 

variance(^^x) = FCt. 

Since C is symmetric, we can write 

max variancef^'^x) = IICIL. 

Using the fact that j|C||2 = cr^ ||(A^^)"^||2 = erg ||A^||2, and Equation ((2Jl) . we get 

max variance (£"^i;) = cr^Krcffo)^ 

IKI|2 = 1 

or, if we call (T(^-'"i;) the standard deviation of P'x, 

max a(i^x) = (JbK,Ls(b)- 

¥\\2 = l 

Note that ab = maxpn^^i cr(^"^e) since V{e) = erg/. 

Remark 2. Matlab proposes a routine LSCOV that computes the quantities y/cil in a vector 
STDX and the mean squared error MSE using the syntax [X,STDX,MSE] = LSCOV(A,B). 
Then the condition numbers Ki{h) can be computed by the matlab expression STDX/sqrt(MSE). 
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3.3. Perturbation on A and b. Wc now provide the expression of the condition number 
provided in Equation (|2.4|1 and in Equation p.5|) in term of statistical quantities. 
Observing the foUowing relations 



fJb \ (m — n) (.V j 




Ci^aleJ{A^A) ^ and cu ^ al \\ej A'^ , 

where Ci is the i-th column of the variance-covariancc matrix, the condition number of Xi given in 
Formula (|2.4p can expressed as 

The quantity will often be estimated by |jr||2 in which case the expression can be simplified 

^1^1 \\C,\\l , c,, 11x11 = 

From Equation (|2.5p . we obtain 
KLs{A,b) 

The quantity cr^ will often be estimated by ||r||2 in which case the expression can be simplified 
KLs{A,b) 

4. Computation with LAPACK. Section [2] provides us with formulas to compute the con- 
dition numbers Ki and kls- As explained in Section [31 those quantities arc intimately interrelated 
with the entries of the variance-covariancc matrix. The goal of this section is to present practical 
methods and codes to compute those quantities efficiently with LAPACK. The assumption made is 
that the LLSP has already been solved with cither the normal equations method or a QR factoriza- 
tion approach. Therefore the solution vector x, the norm of the residual ||r||2, and the R-factor R 
of the QR factorization of A are readily available (we recall that the Cholesky factor of the normal 
equations is the R-factor of the QR factorization up to some signs). In the example codes, we have 
used the LAPACK routine DGELS that solves the LLSP using QR factorization of A. Note that 
it is possible to have a more accurate solution using extra-precise iterative refinement [H]. 

4.1. Variance-covariance computation. We will use the fact that — Ai||2 is an 

unbiased estimate of a^. We wish to compute the following quantities related to the variance- 
covariance matrix C 

• the i-th column d = a'^ef {A^ A)~^ 

• the i-th diagonal element cu = cr^\\e[Am2 

• the whole matrix C 

We note that the quantities Ci, cu, and C are of interest for statisticians. The NAG routine 
F04YAF [12] is indeed an example of tool to compute these three quantities. 
For the two first quantities of interest, we note that 

\\efA^\l^\\R-^e4l and \\ef {A'^ A)-\ ^ \\R-\R~^ e,)\\^ . 
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4.1.1. Computation of the i-th column C'i. Ci can be computed with two n-hy-n trian- 
gular solves 

R^y = e,; and Rz = y. (4.1) 
The i-th column of C can be computed by the following code fragment. 

Code 1: 

CALL DGELS( 'N', M, N, 1, A, LDA, B, LDB, WORK, LWORK, INFO ) 
RESNORM = DNRM2( (M-N), B(N+1), 1) 
SIGMA2 = RESNORM**2/DBLE(M-N) 
E(1:N) = O.DO 
E(I) l.DO 

CALL DTRSV( 'U', 'T', 'N', N-I+1, A(I,I), LDA, E(I), 1) 
CALL DTRSV( 'U', 'N', 'N', N, A, LDA, E, 1) 
CALL DSCAL( N, SIGMA2, E, 1) 

This requires about 2n^ flops (in addition to the cost of solving the linear least squares 
problem using DGELS). 

Cii can be computed by one n-by-n triangular solve and taking the square of the norm of the 
solution which involves about {n — i + 1)^ flops. It is important to note that the larger i, the less 
expensive to obtain cu. In particular ii i = n then only one operation is needed: c„„ = R^n- This 
suggests that a correct ordering of the variables can save some computation. 

4.1.2. Computation of the i-th diagonal element cu. From ca = ct^ ||ef H^, it 
comes that each cu corresponds to the i-th row of R^^. Then the diagonal elements of C can be 
computed by the following code fragment. 

Code 2: 

CALL DGELS( 'N', M, N, 1, A, LDA, B, LDB, WORK, LWORK, INFO ) 
RESNORM = DNRM2((M-N), B(N+1), 1) 
SIGMA2 = RESNORM**2/DBLE(M-N) 
CALL DTRTRI( 'U', 'N', N, A, LDA, INFO) 
DO I=1,N 

CDIAG(I) = DNRM2( N-I+1, A(I,I), LDA) 
CDIAG(I) = SIGMA2 * CDIAG(I)**2 
END DO 

This requires about n'^/S flops (plus the cost of DGELS). 



4.1.3. Computation of the whole matrix C. In order to compute explicity all the 
coefficients of the matrix C, one can use the routine DPOTRI which computes the inverse of a 
matrix from its Cholesky factorization. First the routine computes the inverse of R using DTRTRI 
and then performs the triangular matrix-matrix multiply R~^R~^ by DLAUUM. This requires 
about 2n^/3 flops. We can also compute the variance-covariance matrix without inverting R using 
for instance the algorithm given in [51 p. 119] but the computational cost remains 2n^/3 (plus the 
cost of DGELS). 

We can obtain the upper triangular part of C by the following code fragment. 
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Code 3: 

CALL DGELS( 'N', M, N, 1, A, LDA, B, LDB, WORK, LWORK, INFO ) 
RESNORM = DNRM2((M-N), B(N+1), 1) 
SIGMA2 = RESNORM**2/DBLE(M-N) 
CALL DPOTRI( 'U', N, A, LDA, INFO) 

CALL DLASCL( 'U', 0, 0, N, N, l.DO, SIGMA2, N, N, A, LDA, INFO) 



4.2. Condition numbers computation. For computing Ki{A, b), we need to compute both 
the i-th diagonal element and the norm of the z-th column of the variance-covariance matrix and 
we cannot use direcly Code 1 but the following code fragment 

Code 4: 

ALPHA2 = ALPHA**2 
BETA2 = BETA**2 

CALL DGELS( 'N', M, N, 1, A, LDA, B, LDB, WORK, LWORK, INFO ) 

XNORM = DNRM2(N, B(l), 1) 

RESNORM = DNRM2((M-N), B(N+1), 1) 

CALL DTRSV( 'U', 'T', 'N', N-I+1, A(I,I), LDA, E(I), 1 ) 

ENORM = DNRM2(N, E, 1) 

K = (ENORM**2)*(XNORM**2/ALPHA2+l.dO/BETA2) 
CALL DTRSV( 'U', 'N', 'N', N, A, LDA, E, 1 ) 
ENORM = DNRM2(N, E, 1) 

K = SQRT((ENORM*RESNORM)**2/ALPHA2 + K) 

For computing all the Ki{A,b), we need to compute the columns Ci and the diagonal ele- 
ments Cii using Formula (|3.4p and then we have to compute the whole variance-covariance matrix. 
This can be performed by a slight modification of Code 3. 

When only b is perturbed, then we have to invert R and we can use a modification of Code 2 (see 
numerical example in Section [5. 2p . 

For estimating KLs{^,b), we need to have an estimate of ||^~^||2- The computation of 
[["^"^IL ''^l^ii'^s to compute the minimum singular value of the matrix A (or R). One way is to 
compute the full SVD of A (or R) which requires 0{n^) flops. As an alternative, jj^^^Hj can be 
estimated for instance by considering other matrix norms through the following inequalities 



\\R or \\R ^j|oo can be estimated using Higham modification [UJ p. 293] of Hager's [IH] 
method as it is implemented in LAPACK p] DTRCON routine (sec Code 5). The cost is 0{n'^). 

Code 5: 

CALL DTRCON( 'I', 'U', 'N', N, A, LDA, RCOND, WORK, IWORK, INFO) 
RNORM = DLANTR( 'I', 'U', 'N', N, N, A, LDA, WORK) 
RINVNORM = (l.D0/RNORM)/RCOND 



9 



We can also evaluate \\R ^||, by considering \\R ^||^ since we have 



1-^ ^\\f ~ 11^ ^IIf 



^tr(C), 



where tr(C) denotes the trace of the matrix C, i.e X]"=i '^a- Hence the condition number of the 
least-squares solution can be approximated by 

Then we can estimate Ki5(A, 6) by computing and summing the diagonal elements of C using 
Code 2. 

When only b is perturbed {a = +oo and (3=1), then we get 



KLs[b) ~ . 

This result relates to [21 P- 167] where tr(C) measures the squared effect on the LLSP solution x 
to small changes in b. 

We give in Table 14.11 the LAPACK routines used for computing the condition numbers of 
an LLSP solution or its components and the corresponding number of floating-point operations 
per second. Since the LAPACK routines involved in the covariance and/or LLSP condition 
numbers have their equivalent in the parallel library ScaLAPACK [S], then this table is also 
available when using ScaLAPACK. This enables us to easily compute these quantities for larger 
LLSP. 

Table 4.1 

Computation of least squares conditioning with (Sca)LAPACK 



condition number 


linear algebra operation 


LAPACK routines 


flops count 


b) 


y = Gi and Rz = y 


2 calls to (P)DTRSV 


2n^ 


all 6), i = \,n 


RY = I and compute YY^ 


(P)DPOTRI 


2n^/3 


all Ki{b), i = l,n 


invert R 


(P)DTRTRI 




KLs{A,b) 


estimate j or oo 
compute p 


(P)DTRCON 
(P)DTRTRI 


n'V3 



Remark 3. The cost for computing all the Ki{A,b) or estimating 6) is always 

0{n^). This seems affordable when we compare it to the cost of the least squares solution using 
Householder QR factorization (2mn^ — 2n'^/3) or the normal equations (mn^ -I- n^/3) because we 
have in general m ^ n. 



10 



5. Numerical experiments. 



5.1. Laplace's computation of the mass of Jupiter and assessment of the validity 
of its results. In ^TOj, Laplace computes the mass of Jupiter, Saturn and Uranus and provides 
the variances associated with those variables in order to assess the quality of the results. The 
data comes from the French astronomer Bouvart in the form of the normal equations given in 
Equation (|5.ip . 

f 2602z5 
f 5722z5 

1.3371Z.5 
1.1128Z5 
46.310Z5 
129z5 



7212.600 
-738297.800 
237.782 
-40.335 
-343.455 
-1002.900 

(5.1) 

For computing the mass of Jupiter, we know that Bouvart performed m = 129 observations and 
there are n — 6 variables in the system. The residual of the solution \\b — Ax\\2 is also given by 
Bouvart and is 31096. On the 6 unknowns, Laplace only seeks one, the second variable zi. The 
mass of Jupiter in term of the mass of the Sun is given by zi and the formula: 



795938zo - 12729398zi + 6788. 2z2 - 1959. Oza + 696.13z4 
-12729398zo + 424865729zi - 153106.5z2 - 39749.1^3 - 5459z4 
6788. 2zo - 153106. 5zi + 71.8720z2 - 3. 2252^3 + 1.2484z4 4 
-1959.0zo - 39749. Izi - 3.2252z2 + 57.1911z3 + 3.6213z4 4 
696.1320 - 5459zi 4 1.2484z2 4 3.6213z3 + 21.543z4 4 
2602zo 4 5722zi 4 1.3371z2 + 1. 1128^3 4 46.310z4 



mass of Jupiter 



1 4zi 



1067.09 

It turns out that the first variable zq represents the mass of Uranus through the formula 

1 + 20 



mass of Uranus 



19504' 



If we solve the system ()5.1|) . we obtain the solution vector 
Solution vector 

0.08954 -0.00304 -11.53658 -0.51492 5.19460 -11.18638 
From Zi , we can compute the mass of Jupiter as a fraction of the mass of the Sun and we obtain 

1070. This value is indeed accurate since the correct value according to NASA is 1048. From zq, 

we can compute the mass of Uranus as a fraction of the mass of the Sun and we obtain 17918. 

This value is inaccurate since the correct value according to NASA is 22992. 

Laplace has computed the variance of zq and zi to assess the fact that zi was probably correct 
and Zq probably inaccurate. To compute those variances, Laplace first performed a Cholesky 
factorization from right to left of the system ()5.ip . then, since the variables were correctly ordered 
the number of operations involved in the computation of the variances of zq and zi were minimized. 
The variance-covariance matrix for Laplace's system is: 



/ 0.005245 



V 



-0.000004 -0.499200 0.137212 
0.000004 0.009873 0.003302 
71.466023 -5.441882 
10.860492 



0.235241 
0.002779 
-16.672689 
5.418506 
66.088476 



-0.186069 
-0.001235 
14.922752 
-4.896579 
-28.467391 
15.874809 



Our computation gives us that the variance for the mass of Jupiter is 4.383233 • 10~^. For 
reference, Laplace in 1820 computed 4.383209 • 10~^. (We deduce the variance from Laplace's 
value 5.0778624. To get what we now call the variance, one needs to compute the quantity: 
1/(2 * 10 * *5.0778624) * m/(m - n).) 

From the variance-covariance matrix, one can assess that the computation of the mass of 
Jupiter (second variable) is extremely reliable while the computation of the mass of Uranus (first 
variable) is not. For more details, we recommend to read [18] . 
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5.2. Gravity field computation. A classical example of parameter estimation problem 
is the computation of the Earth's gravity field coefficients. More specifically, we estimate the 
parameters of the gravitational potential that can be expressed in spherical coordinates (r, 9, A) 
by U 

i^a. /^\ ^ + 1 ^ _ _ 

e=0 ^ ^ m=Q 

where G is the gravitational constant, M is the Earth's mass, R is the Earth's reference radius, the 
Pirn represent the fully normalized Legendre functions of degree £ and order m and Cf,ra,Sirn are 
the corresponding normalized harmonic coefficients. The objective here is to compute the harmonic 
coefficients Cf„i and Sim the most accurately as possible. The number of unknown parameters is 
expressed by n = [Imax + 1)^- These coefficients are computed by solving a linear least squares 
problem that may involve millions of observations and tens of thousands of variables. More details 
about the physical problem and the resolution methods can be found in [3] . The data used in the 
following experiments were provided by CNE^ and they correspond to 10 days of observations 
using GRACeH measurements (about 166, 000 observations). We compute the spherical harmonic 
coefficients C^,„ and Sim up to a degree imax = 50; except the coefficients Cn, ^n, Coo, Cio that 
are a priori known. Then we have n = 2, 597 unknowns in the corresponding least squares problems 
(note that the GRACE satellite enables us to compute a gravity field model up to degree 150). The 
problem is solved using the normal equations method and we have the Cholesky decomposition 
A^A ^ U^U. 

We compute the relative condition numbers of each coefficient Xi using the formula 

nt'\b)^\\elU-W\b\\J\x^,l 

and the following code fragment, derived from Code 2, in which the array D contains the normal 
equations A^A and the vector X contains the right-hand side A'^b. 

CALL DPOSV( 'U', N, 1, D, LDD, X, LDX, INFO) 
CALL DTRTRI( 'U', 'N', N, D, LDD, INFO) 
DO I=1,N 

KAPPA(I) = DNRM2( N-I+1, D(I,I), LDD) * BNORM/ABS(X(I)) 
END DO 

Figure 15.11 represents the relative condition numbers of all the n coefficients. We observe 
the disparity between the condition numbers (between 10^ and 10*). To be able to give a 
physical interpretation, we need first to sort the coefficients by degrees and orders as given in the 
development of V{r,9, X) in Expression (|5.2p . 

In Figure [SHI we plot the coefficients C^m as a function of the degrees and orders (the curve with 
the Sim is similar). We notice that for a given order, the condition number increases with the 
degree and that, for a given degree, the variation of the sensitivity with the order is less significant. 
We can also study the effect of regularization on the conditioning. The physicists use in 
general a Kaula [16| regularization technique that consists of adding to A'^ A a diagonal matrix 
D = diag{0, ■ ■ ■ , 0, i5, • • • ,S) where (5 is a constant that is proportional to ^ — and the nonzero 
terms in D correspond to the variables that need to be regularized. An example of the effect 
of Kaula regularization is shown in Figure 15.31 where we consider the coefficients of order also 
called zonal coefficients. We compute here the absolute condition numbers of these coefficients 
using the formula Ki{b) ~ \\^TU^^\\2- ^^te that the Ki{b) are much lower that 1. This is not 
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surprising because typically in our application ||6||2 ~ 10^/ and \xi\ ~ lO""*^^ which would make 
the associated relative condition numbers greater than 1. We observe that the regularization is 
effective on coefficients of highest degree that are in general more sensitive to perturbations. 




500 1000 1500 2000 2500 

# gravity field coefficients 

Fig. 5.1. Amplitude of the relative condition numbers for the gravity field coefficients. 




'^^'^ order 
Fig. 5.2. Conditioning of spherical harmonic coefficients Cim (2<f<50, l<ni< 50). 
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Fig. 5.3. Effect of regularization on zonal coefficients Ciq {2 < £ < 50) 



6. Conclusion. To assess the accuracy of a linear least squares solution, the practitioner of 
numerical linear algebra uses generally quantities like condition numbers or backward errors when 
the statistician is more interested in covariance analysis. In this paper we proposed quantities 
that talk to both communities and that can assess the quality of the solution of a least squares 
problem or one of its component. We provided pratical ways to compute these quantities using 
(Sca)LAPACK and we experimented these computations on pratical examples including a real 
physical application in the area of space geodesy. 
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